
#ifndef __RBFFUNCTION_3TO1_H__
#define __RBFFUNCTION_3TO1_H__

#include "rbfFunction.h"

class RbfFunction3to1 : public RbfFunction {
   
public:
   /// Default constructor
   RbfFunction3to1():
      RbfFunction() {
      npolyterms = 0;
   }
   /// Destructor
   virtual ~RbfFunction3to1() {
   }
   
   void setCentersAndWeights(vector<Point3> &c, vector<Number> &w) {
      centers = c;
      weights = w;
   }
   
   void setEvalPoints(vector<Point3> &points) {
      pPoints = &points;
   }
   
   /// Performs the product: w = iM*fvalues
   /// @param iM the RBF inverse matrix
   /// @param cntrs the centers array
   /// @param w the weights array
   void genWeights(vector<Number> &fvalues) {
      
      assert(centers.size() == fvalues.size());
      unsigned int n = fvalues.size();
      assert(npolyterms == iM.m-n);
      Matrix RightSide(iM.m, 1);
      for (unsigned int i=0; i<n; ++i) {
         RightSide(i,0)=fvalues[i];
      }
      
      /// additional restrictions initialization
      for (unsigned int i=0; i<npolyterms; ++i) {
         RightSide(n+i,0)=0.0;
      }
      Matrix lambdas = iM*RightSide;
      
      weights.resize(iM.m);
      for (unsigned int i=0; i<weights.size(); ++i) {
         weights[i]=lambdas(i,0);
      }
   }
   
   void execute(vector<Number> &values) {
      
      assert(pPoints->size() == values.size());
      for (uint i=0; i<values.size(); ++i) {
         values[i] = evalf((*pPoints)[i]);
      }
   }
   
   Number evalf(const Point3 &p) {
      
      const uint &ncenters = centers.size();
      Number r=0;
      /// radial component
      for (unsigned int i=0; i<ncenters; ++i) {
         const Point3 &P = centers[i];
         const Number &PHI = phi(p-P);
         r += weights[i]*PHI;
      }
      
      /// polynomial component
      /// parte das colunas de lambdas sao os coeficientes do polinomio
      /// que corresponde ao componente x, y ou z
      for (unsigned int i=0; i<npolyterms; ++i) {
         const Number &valpoly = polynomial(p, i);
         r += weights[ncenters+i]*valpoly;
      }
      return r;
   }
   
private:
   vector<Number>    weights;         // the vector of weights
};

#endif
